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ABSTRACT 

Large-scale star formation in disk galaxies is hypothesized to be driven by 
global gravitational instability. The observed gas surface density is commonly 
used to compute the strength of gravitational instability, but according to this 
criterion star formation often appears to occur in gravitationally stable regions. 
One possible reason is that the stellar contribution to the instability has been 
neglected. We have examined the gravitational instability of the Large Magel- 
lanic Cloud (LMC) considering the gas alone, and considering the combination 
of collisional gas and collisionless stars. We compare the gravitationally unsta- 
ble regions with the on-going star formation revealed by Spitzer observations of 
young stellar objects. Although only 62% of the massive young stellar object 
candidates are in regions where the gas alone is unstable, some 85% lie in re- 
gions unstable due to the combination of gas and stars. The combined stability 
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analysis better describes where star formation occurs. In agreement with other 
observations and numerical models, a small fraction of the star formation occurs 
in regions with gravitational stability parameter Q > 1. We further measure the 
dependence of the star formation timescale on the strength of gravitational in- 
stability, and quantitatively compare it to the exponential dependence expected 
from numerical simulations. 

Subject headings: Galaxies: kinematics and dynamics — galaxies: ISM — galax- 
ies: stellar content — Magellanic Clouds — stars: formation 



INTRODUCTION 



Observations of nearby galaxies suggest that there exists a gas surface density thresh- 



old for star formation in a galactic disk (Kennicutt 1989; Martin & Kennicutt 2001). This 



threshold is hypothesized to be formed by the threshold for gravitational stability of the disk, 



although alternative explanations have been offered (e.g., Schaye 2004). Gravitationally un- 



stable regions are commonly identified by the Toomre criterion (Toomre 1964 Goldreich 



& Lynden-Bell 1965), which considers the stability of a single-component disk against ax- 



isymmetric perturbations. However, this criterion, when applied to observed gas surface 
densities, does not always correctly predict the locations of star formation (e.g., Martin & 



Kennicutt 2001 Wong & Blitz 2002). Stellar contribution to the gravitational instability 



could be important; for example, the giant molecular clouds in the Large Magellanic Cloud 
(LMC) have been shown to correlate with stellar surface density, indicating the increase of 



star formation towards regions of higher stellar gravity (Fukui 2007). 



To crudely include the stellar contribution, Jog & Solomon (1984) derived the stability 



criterion for a two-fluid disk. This stability criterion differs from the Toomre criterion by 
only a linear correction factor, when the gas-to-stellar ratios of surface density and velocity 
dispersion are constant across a galaxy, or when the gas and stellar components are only 



weakly coupled (Wang & Silk 1994). To correctly include the stellar contribution, Gammie 



(1992) and Rafikov (2001 ) analyzed the stability of a composite disk consisting of a collisional 



gas component and a population of collisionless stellar components. The resulting instability 



criterion has been used to explain the dust morphology of edge-on disk galaxies (Dalcanton, 



Yoachim, & Bernstein 2004). 



We study the relationship between gravitational stability and star formation in the 
LMC. The LMC is chosen because high linear resolution can be obtained (1" corresponds 



to 0.25 pc at 50 kpc, Feast 1999), its internal and foreground extinctions are small, and 
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its moderate inclination allows spatial and kinematical mapping of the disk with minimum 
confusion along the line of sight. Both stars and the interstellar medium (ISM) of the 
LMC have been extensively surveyed. Recent Spitzer Space Telescope observations have 
been used to identify massive young stellar objects (YSOs), which mark the sites of current 
star formation (R. A. Gruendl et al., in preparation; B. A. Whitney et al., in preparation). 
Massive YSOs provide better probes than H II regions to study the relationship between the 
ISM and star formation because their stellar energy feedback has not significantly altered 
the distributions and physical conditions of the ambient ISM. 

In this paper, we investigate the relationship between global star formation and grav- 
itational instability in the LMC using two different stability analyses: assuming that the 
gaseous disk is decoupled from the stellar disk and considering only the gravity of the gas 



2.1 ), or assuming a collisionless stellar disk and a collisional gas disk and using the Rafikov 



(2001) stability criterion (§ 2.2). The results are discussed in §w 



GRAVITATIONAL INSTABILITY 



2.1. Gas Alone 



A thin, differentially rotating gaseous disk is unstable against axisymmetric perturba- 
tions when 



Qg = 



KCn 



< 1, 



where Q g is generally referred to as the Toomre Q parameter for a gaseous disk, k is the 
epicyclic frequency, c g is the isothermal sound speed of the gas, G is the gravitational con- 
stant, and E 9 is the unperturbed surface density of the gas ( Goldreich fc Lynden-Bell||1965 ). 



To apply this analysis to a disk galaxy, approximations need to be made. First, rather 
than consisting of a single isothermal component, the real ISM consists of gas with temper- 
atures ranging from 10 K to 10 7 K in multiple phases including cold atomic and molecular, 
warm atomic and ionized, and hot ionized gas. However, the gas mass is dominated by cold 
atomic and molecular material in turbulent motion, except possibly in star-forming regions 
where much of the ISM is ionized. Therefore, we consider only the neutral atomic and 
molecular gas mass. The atomic medium of galaxies generally has nearly constant velocity 



dispersion (van der Kruit & Shostak 1982 Shostak & van der Kruit 1984 Dickey, Hanson, 



& Helou 1990 Petric & Rupen 2007), and has a vertical thickness that is well modeled by 



assuming hydrostatic equilibrium with the gravitational potential at the velocity dispersion 
(e.g., Malhotra 1995). We assume the gas disk of the LMC has an effective sound speed 



representative of the velocity dispersion c„ = 5 km s 1 (Dib, Bell, & Burkert 2006 Petric 



& Rupen 2007). Second, £ 9 in equation ([I]) is the unperturbed gas surface density, while 



the observed gas surface density is perturbed. The surface density of a gravitationally un- 
stable region grows exponentially as it collapses; thus, using the observed surface density in 
equation ([TJ will correctly diagnose gravitationally unstable regions, but Q g may be under- 
estimated if collapse has set in. Therefore, we adopt the observed gas surface density for E s 
in equation Q to identify gravitationally unstable regions. 

The neutral atomic and molecular components of the ISM in the LMC have been well 
surveyed. For the neutral atomic component, we use the H I column density map derived 
from the combined data set of the Australia Telescope Compact Array and the Parkes 
multibeam receiver (Kim et al. 2003). This H I map has a resolution of ~1 , or 15 pc. 



For the molecular component, we use the NANTEN CO survey of the LMC conducted by 
Fukui et al. (1999, 2001). These maps have a resolution of ~2'6, or 40 pc. We scale maps 
of integrated column density by the CO-to-H2 conversion factor, X = 5.4 x 10 20 H2 atoms 
~ 2 (K km s -1 ) -1 (Blitz et al. 2007) to obtain the H2 mass distribution. Both the H I and 



cm 



H 2 maps are clipped at the la level to prevent negative densities and minimize the effects 
of noise. Furthermore, we account for He and heavier elements by assuming that the mass 
fraction of H is X ~ 0.7. From the H I map, we find a total mass for the neutral atomic 



component of 5.1 x 10 8 M , consistent with the value found by Kim et al. (1998). For the 
molecular component, we find a total mass of 6.6 x 10 7 M , which is on the high side of the 



4-7xl0 7 M Q given by Mizuno et al. (1999, 2001), as we used a higher CO-to-H 2 conversion 



factor and included He and heavy elements. We add these two maps to obtain a map of 
total gas surface density with a resolution of 40 pc pixel -1 . The resulting map is shown in 
Figure [jji. 

The remaining parameter needed in equation (pj) is the epicyclic frequency k, defined 

by 



K 2 = 



2V 2 



VdR 



where R is the galactocentric distance and V = V(R) is the circular velocity as a function of 



R. To calculate k at a given R in the LMC, we use the rotation curve derived by Kim et al. 



(1998). The rotation curve is a best-fit to the H I and the carbon star (Kunkel et al. 1997) 



measurements inside and outside of about 3.2 kpc, respectively. We find the derivatives 
of the circular velocity and thus the epicyclic frequency by central finite differences with 
an increment of 250 pc, then we use a cubic spline to interpolate the discrete results to a 
continuous domain. Finally, the center and the orientation of the disk of the LMC need 
to be specified. To be consistent with the H I kinematics, we adopt the same center and 



orientation determined by Kim et al. (1998). The kinematic center of the LMC disk is at 
a = 05 h 17 m 6 and 5 = -69°02' (J2000) with an inclination of i = 33° and a line of nodes 
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at a P.A. of —12°. Using this geometry, the total gas surface density map is deprojected to 
determine the galactocentric distance of each pixel, allowing us to construct a map of Q g 
calculated from equation (fi]). This map is shown in Figure [T}d. 

The contours in Figure [TJd delineate the critical boundary Q g = 1 that surrounds regions 
gravitationally unstable due to gas alone. Also marked in the figure are candidate massive 
YSOs from R. A. Gruendl et al. (in preparation), indicating regions of on-going star forma- 
tion. In the top panel of Figure |2j we show the number distribution of YSO candidates and 
the corresponding cumulative fraction with respect to Q g . While ~62% (153 of 245) of the 
YSO candidates do lie in regions with Q g < 1, the rest are distributed in regions with values 
as high as Q g ~ 3. 

Some of the candidates with Q g < 1 might be argued not to result from large-scale 
gravitational instabilities because the unstable regions in which they reside are too small. If 
we exclude regions < 100 pc in diameter, only 47% (116) of the YSO candidates would be 
associated with gravitationally unstable regions by this criterion. A frequently used empirical 



correction for the contributions from the stellar component (Jog & Solomon 1984 Wang & 



Silk|1994 ) and the effect of the disk scale height ( Jog fc Solomon|[l984 ) is to raise the critical 



value to Q g ~ 1.4 (Kennicutt 1989; Martin & Kennicutt 2001). Even if we apply this linear 
correction factor to the Toomre criterion, ~25 % of the YSO candidates still appear to 
reside in stable regions. Considering the stability of only the gas disk or applying a simple 
linear correction to the Toomre criterion may not be sufficient to account for star formation 
activity. 



2.2. Gas and Stars Together 



We now also consider the contribution of stars to gravitational instability. To include 
them, we follow Rafikov's (2001) treatment of a disk galaxy consisting of a collisional gas 
disk and a collisionless stellar disk. The instability condition becomes 

12 1 

Q s q 



Q 



1 - e~*I a (q 2 ) 



q 



> 1. 



(3) 



In equation ([3]), Q g is the stability parameter for the gas derived by Goldreich & Lynden-Bell 
(1965), as defined in equation ([!]), while 



Qs 



(4) 



is the stability parameter for the stars derived by Toomre (1964), where a s is the stellar 
radial velocity dispersion, S s is the stellar surface density; R = c g /a s and q = ka s j k are two 
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dimensionless parameters, with k being the wavenumber of the axisymmetric perturbations; 
and Jq is the Bessel function of order zero. 



Similarly to our procedure in § 2.1 we use local observed values to evaluate Q sg . To 



estimate the stellar surface density E s , we use the number density of red giant branch (RGB) 
and asymptotic giant branch (AGB) stars, as they are part of the old stellar population and 
therefore trace the overall mass distribution of the stellar disk. RGB and AGB stars are 
luminous and distinctive in the [J — K s ] vs. K s color-magnitude diagram and much less 
confused by faint background galaxies or Galactic stars. To select these stars, we follow a 



procedure similar to that outlined by van der Marel (2001) but use only the Two Micron 
All Sky Survey Point Source Catalogue (Skrutskie et al. 2006) and the criteria: [J — K s ] < 



(22— [if s ]) /10.5 and [K s ] < 14.5. The source counts are then binned with a resolution of 40 pc 
pixel -1 and Gaussian smoothed with a dispersion of 100 pc. The smoothing scale is chosen 
to minimize pixel-to-pixel variations due to the small number density of the tracer RGB and 
AGB stars, while retaining good spatial resolution for the analysis; this scale is also small 
compared to the large-scale gravitational instability in the analysis. In the resulting map, 
a faint Galactic background with a gradient is present. Therefore, we fit and subtract this 
background with exponential and linear dependencies on the Galactic latitude and longitude, 
respectively. Finally, by adopting the total stellar mass of 2 x 10 9 M & estimated by Kim et 



al. ( 1998[ ), we normalize this map to obtain the absolute stellar surface density map shown 
in Figure [1J (cf. Fig. 2 of van der Marel 2001). We adopt a constant stellar radial velocity 
dispersion of 15 km s _1 , after considering the velocity dispersions of carbon stars (Kunkel et 



al. 1997 Alves k Nelson 2000 Graff et al. 2000; van der Marel et al. 2002), red supergiants 



(Prevot et al. 1989 Olsen & Massey 2007), and young globular clusters (Freeman et al. 



1983). 



With these values assessed, the parameter Q sg defined in equation ^ now depends only 
on the perturbation wavenumber k, or equivalently, the perturbation wavelength A = 2ir/k 
at each pixel (see Fig. [3J. Since Q sg decreases with the degree of gravitational instability 
and our objective is to locate gravitationally unstable regions, we define its value at each 
pixel as its global minimum with respect to A. The derived Q sg map of the LMC is shown 
in Figure [T]i, with the identified YSOs overlaid. 

Each point is unstable to perturbations over a finite range of wavelengths, as shown 
by Figure [3} If the size of a gravitationally unstable region is smaller than the minimum 
unstable wavelength A min , no perturbation inside the region can have long enough wavelength 
to drive an instability. Therefore, we also find A m i n in each pixel, and compare this value 
with the size of each region containing YSO candidates by visual inspection. The contours 
of Q sg = 1 are plotted over grayscale presentations of A m j n in Figure |4| so that the sizes of 
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the Q sg < 1 regions can be compared with the minimum unstable wavelengths directly. 

Figure [TJd shows that, because of the stellar contribution, a much larger fraction of the 
galaxy is unstable than would have been deduced just from the gas contribution shown in 
Figure [T]d. The vast majority of the massive YSO candidates are in fact located within 



gravitationally unstable regions, in sharp contrast to the results in § 2.1 where only the 
gas is considered. Particularly noticeable is the vicinity of the LMC bar region, where on- 
going star formation is observed. This bar region, having low gas surface densities and high 
epicyclic frequencies, is stable when only the gas disk is considered, but becomes unstable 
when the stellar component is added. The bottom panel of Figure [2] shows the number 
distribution of YSO candidates and the corresponding cumulative fraction with respect to 
Q sg . About 86% (212 of 245) of the YSO candidates are located in regions where Q sg < 1. 
Figure [4] shows 12 of the 212 YSOs in probably stable regions having Q sg < 1 but region size 
< A m i n . Removing these 12 still leaves 82% of the YSOs whose formation is attributed to 
gravitational instability of the disk of stars and gas. One caution is that finite thickness of 
the disk may lower the critical Q sg value slightly, depending on the scale height of the disk 
relative to the unstable wavelengths ( Jog fc Solomon|[l984 ). Given an estimated scale height 
of ~180 pc (Kim et al. 1999) and kiloparsec-scale unstable wavelengths in the LMC, this 
effect may be small. The results demonstrate that the stellar contribution to gravitational 
instability can be significant. Gravitational instability of the full disk appears responsible 
for most large-scale star formation activity. 



DISCUSSION AND CONCLUSIONS 



Schaye (2004) has argued that a constant gas surface density threshold is a good in- 
dicator in predicting the edge of the star forming disk. He suggested that, for the LMC 
disk, this threshold is ~4 M pc -2 . Figure [T^, shows the contours of 4 M pc -2 in total 
neutral gas surface density. The contours encompass the majority of the LMC disk except 
the central cavities of the superg iant shells LMC-4 and LMC-8 ( |Meaburn|[l980l ) . Clearly, all 



the massive YSO candidates are located where the gas surface density is greater than 4 M 
pc~ 2 , but this constant surface density threshold does not account for the distribution of 
the star forming sites, especially in the inner parts of the LMC disk. Therefore, Schaye's 
(2004) constant surface density threshold is a necessary but not sufficient condition for star 
formation. Detailed analysis of the large-scale gravitational instability, as shown in § [2j is 
still needed. 

Another measure of the response of star formation to disk instability is the relation 
between Q sg and the effective timescale for star formation r s f. This can be approximated if 
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we assume the average lifetime of massive YSOs is constant. We normalize the number of 
YSOs, N*, within a region having a given range of Q sg by the clX6cl ? clS represented by the 
total number of pixels iVp ix within that region. Figure [5] presents N*/N P i x , which is oc t~^, 
for a number of Q sg bins. The errors in the ratio are estimated assuming Poisson statistics 
on the number of YSOs and number of pixels, while the horizontal error bars show the 
size of the bins. A clear log-linear relationship is found, implying that the response of star 
formation depends exponentially on the value of Q sg . Li, Mac Low, & Klessen (2005) also 



found an exponential dependence r s / cx exp(aQ sg ), with a = 4.2 ±0.3 from global numerical 
simulations of galactic disks. A linear fit to the data in Figure [5] gives a slope that can be 
interpreted as a local value of a = 2.7 ± 0.2. This behavior appears to be a general property 
of nonlinear gravitational instability. 

We have made the simple assumption that the structure of the LMC is a thin, flat disk 
with negligible stellar halo. In reality, tidal features due to galaxy-galaxy interactions do 



exist in the outer parts of the disk (e.g., van der Marel 2001 ). Furthermore, at some locations, 
there exist multiple H I velocity components, which we have simply summed to obtain total 
gas surface density. Other uncertainties stem from our assumed constant gas effective sound 
speed (or velocity dispersion) and stellar radial velocity dispersion. These affect the relative 
contribution of each component to the gravitational instability: the higher the local velocity 
dispersion, the less contribution a component has. The sound speed of the gas we adopt, 
5 km s -1 , is a typical average for the neutral gas. The radial velocity dispersion of stars is 
also uncertain. As gravitational instability is much more sensitive to the gas component due 
to its low sound speed, it is less affected by the uncertainty in the stellar velocity dispersion. 
Yet another effect to consider is the scale height of the disk, which tends to stabilize the 



disk (Toomre 1964 Jog & Solomon 1984). Given all these uncertainties, nevertheless, the 
exponential dependence of star formation rate on Q sg that we find appears to be robust. 

From our analysis of the gravitational instability of the LMC, we conclude that the 



contribution of the stellar disk to gravitational instability cannot be ignored (Jog & Solomon 



1984 


Gammie 


1992 


Rafikov 


2001 



that ~85% of the massive YSO candidates lie in gravitationally unstable regions, implying 
that star formation occurs predominantly in these regions. It appears that gravitational 



instability of the disk drives most large-scale star formation, as proposed by Elmegreen 



(2002 


), 


Kravtsov 


(2003 


), and 


Li et al. 


(2005 
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Fig. 1. — (a) Total gas surface density distribution in the Large Magellanic Cloud. Only 
neutral species are included in this map. The grayscale runs linearly from to 100 M pc~ 2 . 
The contour lines show a surface density of 4 M & pc~ 2 . (b) Comparison between massive 
star-forming sites traced by young stellar object candidates (red dots) and regions where 
the gas alone is gravitationally unstable. The Toomre parameter for the gas Q g is shown in 
grayscale, running logarithmically from 5.0 to 0.2. The solid lines delineate the critical value 
of Q g = 1 inside which the regions are gravitationally unstable, (c) Total stellar surface 
density distribution. The grayscale runs linearly from to 200 M pc~ 2 . (d) As in (b), but 
for regions where the gas and stars together are gravitationally unstable, rather than just 
the gas alone. 
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Fig. 2. — Number distribution of YSO candidates with respect to the Toomre stability pa- 
rameter of the pixel in which they are located. The top and bottom panels show respectively 
the Toomre parameters for the gas alone Q g , and the stars and gas together Q sg . The solid 
lines show the cumulative fractions of YSOs with decreasing Q. The vertical dotted lines 
denote the critical value Q = 1. 
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Perturbation Wavelength X 



Fig. 3. — Schematic plot of the Q sg parameter as a function of wavelength A of radial 
perturbations. The disk of stars and gas is stable against perturbation of a given A when 
Q sg > 1; otherwise, it is unstable. We measure the value of Q sg at its minimum in each 
pixel, occurring at A = A crit (point A). We also find the minimum wavelength A min such that 
the disk is marginally unstable (point B). 
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Fig. 4. — Map of the minimum unstable wavelength A m ; n , shown in grayscale. The solid 
contours delineate the critical boundary of Q sg = 1, while the red dots mark the massive 
YSO candidates. A m i n is defined by A cr ;t (see Fig. [3]) when Q sg > 1. 
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Fig. 5. — Ratio of the number of YSO candidates iV* to the number of pixels iVp ix in each 
Q sg bin, denoted by squares. The vertical error bars are estimated by Poisson statistics, 
while the horizontal error bars show the bin sizes. The solid line is the best-fit, of slope 
-2.7 ±0.2. 



